function delpf = Newton2wtV_2(model,fid,kbasis,cgp)

% NEWTON2 Calculates minimizer of line search with quadratic approx
%
% Expects MODEL.{PVEC,PVECTMP}, where PVECTMP represents the direction
% of search.  Also expects FID.{SIG,EST}, KBASIS, and CGP reconstruction
% parameter structure.

delta = 1e-6;

if isfield(fid,'N') 
    NF = fid.N; %this one
else
    NF = length(fid.sig1); 
end
if isfield(model,'N') 
%     NM = model.N; 
else
%     NM = length(model.pvec)/3; %and this one
end
perturb = delta*model.pvectmp;  % from DJ
modelDelta = rmfield(model,'pvectmp');
kmsk1=(fid.estwt);%windowing for signal model
% no perturbation
fid.sig1 = fid.sig1(1:NF);
fid.est1 = fid.est1(1:NF);
fid.sig2 = fid.sig2(1:NF);
fid.est2 = fid.est2(1:NF);
f_0 = jeval_2(fid,modelDelta,cgp);

% perturb one way
modelDelta.pvec = model.pvec + perturb;
fid.est1 = kmsk1(1:fid.N).*fid1syn(modelDelta,kbasis,NF);
fid.est2 = kmsk1(1:fid.N).*fid2syn(modelDelta,kbasis,NF);
f_a_p = jeval_2(fid,modelDelta,cgp);

% perturb the other way
modelDelta.pvec = model.pvec - perturb;
fid.est1 = kmsk1(1:fid.N).*fid1syn(modelDelta,kbasis,NF);
fid.est2 = kmsk1(1:fid.N).*fid2syn(modelDelta,kbasis,NF);
f_a_n = jeval_2(fid,modelDelta,cgp);

f_1_0 = (f_a_p-f_a_n)/(2*delta); % First derivative @ alpha = 0
f_1_a_p = (f_a_p-f_0)/delta; % First derivative @ alpha = delta/2
f_1_a_n = (f_0-f_a_n)/delta; % First derivative @ alpha = -delta/2
f_2_0 = (f_1_a_p-f_1_a_n)/delta; % Second derivative @ alpha = 0

delpf1 = max(0,-f_1_0/f_2_0); % min point

delpf = delpf1;

% an experiment to determine whether further refinement might help.  The
% conclusion is that the quadratic approximation above nails it pretty
% well.  The typical delpf gets smaller with iterations, so delta is too
% small at first but a little too big later.  However, the minimum is very
% close to the refined minimum regardless.
%
% if ~(delpf1 == 0)
%     % perturb to estimated minimum
%     modelDelta.pvec = model.pvec + delpf1*model.pvectmp;
%     fid.est = fidsyn(modelDelta,kbasis,NF);
%     f_a_del = jeval(fid,modelDelta,cgp);
% 
%     % perturb double
%     modelDelta.pvec = model.pvec + 2*delpf1*model.pvectmp;
%     fid.est = fidsyn(modelDelta,kbasis,NF);
%     f_a_2del = jeval(fid,modelDelta,cgp);
% 
%     delpf = delpf1*(-4*f_a_del + 3*f_0 + f_a_2del)/(-2*f_a_del + f_0 + f_a_2del)/2;
% %     fm = min([f_a_n f_0 f_a_p f_a_del f_a_2del]);
% % plot([-delta 0 delta],[f_a_n f_0 f_a_p]-fm,'o',[0 delpf1 2*delpf1],[f_0 f_a_del f_a_2del]-fm,'s',delpf,0,'x',delpf1,0,'d')
% % drawnow
% % pause(1)
% end
